Se você gostou do Kernel não esqueça de deixar um upvote!

Carregar dependência:

knitr::opts_chunk$set(echo = TRUE, fig.width = 11, warning = F, error = F, message = F)

library(tidyverse)
library(tidymodels)
library(vip)

# número de processadores para paralelizar
ncores <- 4

1 Definindo o problema

Este kernel tem como objetivo praticas o ajuste do modelo xgboost: eXtreme Gradient Boosting utilizando o framework tidymodels.



Fonte: https://xgboost.ai/

Fonte: https://www.tidymodels.org/

2 Leitura dos dados

2.1 Informações sobre os dados

Listagem de atributos (traduzido livremente):

feature descrição
idade contínua
classe de trabalho Particular, Auto-emp-não-inc, Auto-emp-inc, Federal-gov, Local-gov, State-gov, Sem salário, Nunca trabalhou
fnlwgt contínuo
educação Bacharelado, Ensino Médio, 11º, Ensino Médio, Prof-escola, Assoc-acdm, Assoc-voc, 9º, 7º-8º, 12º, Mestrado, 1º-4º, 10º, Doutorado, 5º-6º, Pré-Escola
número da educação contínuo
Estado civil cônjuge civil, divorciado, nunca casado, separado, viúvo, cônjuge ausente, cônjuge afetivo
ocupação Suporte técnico, Reparação de ofícios, Outros serviços, Vendas, Gerenciamento executivo, Especialidade prof, Limpadores de manipuladores, Inspeção de máquinas, Adm-administrativo, Agricultura, Pesca, Movimentação de transportes, Casa particular serviço, serviço protetor, forças armadas
Relacionamento Esposa, Filho próprio, Marido, “Não família”, Outro parente, Solteiro
raça branco, asiático-pac-islander, amer-indian-esquimó, outros, preto
sexo Feminino, Masculino
ganho de capital contínuo
perda de capital contínua
horas por semana contínua
país de origem Estados Unidos, Camboja, Inglaterra, Porto Rico, Canadá, Alemanha, Estados Unidos (Guam-USVI-etc), Índia, Japão, Grécia, Sul, China, Cuba, Irã, Honduras, Filipinas, Itália , Polônia, Jamaica, Vietnã, México, Portugal, Irlanda, França, República Dominicana, Laos, Equador, Taiwan, Haiti, Colômbia, Hungria, Guatemala, Nicarágua, Escócia, Tailândia, Iugoslávia, El-Salvador, Trinadad e Tobago,

2.2 Leitura dos dados disponíveis

Leitura conjunto de dados que serão utilizados para treinar o modelo:

adult <- readRDS("input/adult.rds")
adult
## # A tibble: 32,561 x 16
##      age workclass fnlwgt education education_num marital_status occupation
##    <dbl> <chr>      <dbl> <chr>             <dbl> <chr>          <chr>     
##  1    39 State-gov  77516 Bachelors            13 Never-married  Adm-cleri…
##  2    50 Self-emp…  83311 Bachelors            13 Married-civ-s… Exec-mana…
##  3    38 Private   215646 HS-grad               9 Divorced       Handlers-…
##  4    53 Private   234721 11th                  7 Married-civ-s… Handlers-…
##  5    28 Private   338409 Bachelors            13 Married-civ-s… Prof-spec…
##  6    37 Private   284582 Masters              14 Married-civ-s… Exec-mana…
##  7    49 Private   160187 9th                   5 Married-spous… Other-ser…
##  8    52 Self-emp… 209642 HS-grad               9 Married-civ-s… Exec-mana…
##  9    31 Private    45781 Masters              14 Never-married  Prof-spec…
## 10    42 Private   159449 Bachelors            13 Married-civ-s… Exec-mana…
## # … with 32,551 more rows, and 9 more variables: relationship <chr>,
## #   race <chr>, sex <chr>, capital_gain <dbl>, capital_loss <dbl>,
## #   hours_per_week <dbl>, native_country <chr>, resposta <chr>, id <int>

A base possui 32561 linhas e 16 colunas.

3 Framework tidymodels

Aplicaremos a seguir o framework tidymodels

3.1 Amostragem

Separar os dados em treino e teste:

set.seed(32)
adult_split <- initial_split(adult, prop = 0.8, strata = resposta)

3.2 Pré-processamento

Todo o tratamento dos dados será armazenado em uma “receita” de forma que este objeto possa ser executado tanto na hora do treinamento do modelo quanto na hora de testar para submissão.

Obs.: Todos os insights para transformações nos dados foram avaliados após a inspeção do relatório automático desenvolvido utilizando a função DataExplorer::create_report() do pacote DataExplorer

 adult_recipe <- 
  training(adult_split) %>% 
  recipe(resposta ~ .) %>% 
  step_rm(id, education) %>% 
  step_mutate(
    native_country = case_when(
      native_country == "United-States" ~ "USA",
      TRUE ~ "other"),
    occupation = case_when(
      is.na(occupation) ~ "Unemployed",
      TRUE ~ as.character(occupation)),
    workclass = case_when(
      is.na(workclass) ~ "Unemployed",
      TRUE ~ as.character(workclass)),
    capital = capital_gain + capital_loss
  ) %>% 
  step_rm(capital_gain, capital_loss)%>% 
  step_string2factor(all_nominal()) %>%
  step_log(fnlwgt, age) %>%
  step_normalize(fnlwgt, age) %>%
  step_zv(all_predictors()) %>%
  step_center(all_numeric()) %>%
  step_scale(all_numeric()) %>%
  step_dummy(all_nominal(), -all_outcomes())

# bake(prep(adult_recipe), training(adult_split))

Após definir todo o tratamento podemos dar inicio ao workflow:

adult_workflow <- 
  workflow() %>% 
  add_recipe(adult_recipe)

3.3 Validação Cruzada

Agora que nossos dados já foram divididos em treino e teste e já temos a “receita”, vamos especificar como será a validação cruzada:

set.seed(123)
adult_vfold <- vfold_cv(training(adult_split), v = 5, strata = resposta)
adult_vfold
## #  5-fold cross-validation using stratification 
## # A tibble: 5 x 2
##   splits               id   
##   <named list>         <chr>
## 1 <split [20.8K/5.2K]> Fold1
## 2 <split [20.8K/5.2K]> Fold2
## 3 <split [20.8K/5.2K]> Fold3
## 4 <split [20.8K/5.2K]> Fold4
## 5 <split [20.8K/5.2K]> Fold5

3.4 Especificações dos Modelos

Os dois modelos que serão ajustados:

  • rpart: Ajuste de uma árvore (baseline)
  • xgboost: Combinação do ajuste de várias árvores

3.4.1 Baseline - Árvore de Decisão

Definindo o modelo:

adult_rpart_model <- 
  decision_tree(
    min_n = tune(),
    cost_complexity = tune(), 
    tree_depth = tune()) %>%
  set_mode("classification") %>%
  set_engine("rpart")
adult_rpart_model
## Decision Tree Model Specification (classification)
## 
## Main Arguments:
##   cost_complexity = tune()
##   tree_depth = tune()
##   min_n = tune()
## 
## Computational engine: rpart

Definir parâmetros para tuning:

rpart_params <- parameters(
  min_n(),
  cost_complexity(), 
  tree_depth()
)
rpart_params
## Collection of 3 parameters for tuning
## 
##               id  parameter type object class
##            min_n           min_n    nparam[+]
##  cost_complexity cost_complexity    nparam[+]
##       tree_depth      tree_depth    nparam[+]

Montar grade de pesquisa

Alguns métodos para pesquisa em grade são fornecidos pelo pacote dials, para ajustar este modelo baseline utilizaremos o grid_max_entropy():

set.seed(123)
rpart_grid <- grid_max_entropy(rpart_params, size = 20) # local=5 , kaggle=20

Uma representação gráfica da grade:

plotly::plot_ly(rpart_grid, x = ~min_n, y = ~cost_complexity, z = ~tree_depth,
                type = "scatter3d", mode = "markers")

Preparar workflow para rpart:

workflow_adult_rpart_model <- 
  adult_workflow %>% 
  add_model(adult_rpart_model)

Tuning do modelo:

doParallel::registerDoParallel(ncores)

rpart_tune <- 
  workflow_adult_rpart_model %>% 
  tune_grid(
    resamples = adult_vfold,
    grid = rpart_grid,
    control = control_grid(save_pred = TRUE, verbose = F, allow_par = T),
    metrics = metric_set(roc_auc, accuracy, precision, recall)
  )

doParallel::stopImplicitCluster()

Como cada métrica se saiu nos ajustes:

collect_metrics(rpart_tune) %>%
  filter(.metric == "roc_auc") %>%
  select(mean, cost_complexity:min_n) %>%
  pivot_longer(cost_complexity:min_n,
               values_to = "value",
               names_to = "parameter"
  ) %>%
  ggplot(aes(value, mean, color = parameter)) +
  geom_point(alpha = 0.8, show.legend = FALSE) +
  facet_wrap(~parameter, scales = "free_x") +
  labs(x = NULL, y = "AUC")

Melhor combinação:

show_best(rpart_tune, "roc_auc") # top 5
## # A tibble: 5 x 8
##   cost_complexity tree_depth min_n .metric .estimator  mean     n std_err
##             <dbl>      <int> <int> <chr>   <chr>      <dbl> <int>   <dbl>
## 1  0.000000000266         11    18 roc_auc binary     0.897     5 0.00188
## 2  0.00000192             15    25 roc_auc binary     0.896     5 0.00281
## 3  0.00000000708           8    30 roc_auc binary     0.896     5 0.00136
## 4  0.000000000178         14    36 roc_auc binary     0.895     5 0.00259
## 5  0.00000439             12    15 roc_auc binary     0.894     5 0.00201
rpart_best_auc <- select_best(rpart_tune, "roc_auc") # melhor de todos
rpart_best_auc
## # A tibble: 1 x 3
##   cost_complexity tree_depth min_n
##             <dbl>      <int> <int>
## 1  0.000000000266         11    18

Finalizar workflow:

rpart_workflow_final <- finalize_workflow(
  workflow_adult_rpart_model,
  rpart_best_auc
)
rpart_workflow_final
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: decision_tree()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 10 Recipe Steps
## 
## ● step_rm()
## ● step_mutate()
## ● step_rm()
## ● step_string2factor()
## ● step_log()
## ● step_normalize()
## ● step_zv()
## ● step_center()
## ● step_scale()
## ● step_dummy()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Decision Tree Model Specification (classification)
## 
## Main Arguments:
##   cost_complexity = 2.65603871411766e-10
##   tree_depth = 11
##   min_n = 18
## 
## Computational engine: rpart

Atributos mais importantes no ajuste do modelo:

rpart_workflow_final %>%
  fit(training(adult_split)) %>%
  pull_workflow_fit() %>%
  vip(geom = "col")

Finalizar modelo:

rpart_final <- last_fit(rpart_workflow_final, adult_split)
collect_metrics(rpart_final)
## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.846
## 2 roc_auc  binary         0.898

3.4.2 Final Model - XGBoost

Definir modelo

adult_xgb_model <- 
  boost_tree(
  trees = tune(), learn_rate = tune(), # n arvores e cautela para incluir novas arvores
  tree_depth = tune(), min_n = tune(), # parametros da arvore
  loss_reduction = tune(), # gama = 0 arvore complexa, gama = 1 arvore cotoco
  sample_size = tune(), mtry = tune(), # sorteio linhas/colunas
  ) %>% 
  set_mode("classification") %>% 
  set_engine("xgboost",
              lambda = 0,
             nthread = ncores) # n de nucleos 
adult_xgb_model
## Boosted Tree Model Specification (classification)
## 
## Main Arguments:
##   mtry = tune()
##   trees = tune()
##   min_n = tune()
##   tree_depth = tune()
##   learn_rate = tune()
##   loss_reduction = tune()
##   sample_size = tune()
## 
## Engine-Specific Arguments:
##   lambda = 0
##   nthread = ncores
## 
## Computational engine: xgboost

Definir parâmetros para tuning:

xgboost_params <- parameters(
  trees(), learn_rate(),
  tree_depth(), min_n(), 
  loss_reduction(),
  sample_size = sample_prop(), finalize(mtry(), training(adult_split))  
)
# diminuir tempo de processamento local:
xgboost_params <- xgboost_params %>% update(trees = trees(c(100, 500))) 
# tunar lambda:
# xgboost_params <- xgboost_params %>% update(lambda = dials::penalty())
                                            
xgboost_params
## Collection of 7 parameters for tuning
## 
##              id parameter type object class
##           trees          trees    nparam[+]
##      learn_rate     learn_rate    nparam[+]
##      tree_depth     tree_depth    nparam[+]
##           min_n          min_n    nparam[+]
##  loss_reduction loss_reduction    nparam[+]
##     sample_size    sample_size    nparam[+]
##            mtry           mtry    nparam[+]

Preparar workflow para xgboost

workflow_adult_xgb_model <- 
  workflow() %>% 
  add_model(adult_xgb_model) %>% 
  add_recipe(adult_recipe)
workflow_adult_xgb_model
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: boost_tree()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 10 Recipe Steps
## 
## ● step_rm()
## ● step_mutate()
## ● step_rm()
## ● step_string2factor()
## ● step_log()
## ● step_normalize()
## ● step_zv()
## ● step_center()
## ● step_scale()
## ● step_dummy()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Boosted Tree Model Specification (classification)
## 
## Main Arguments:
##   mtry = tune()
##   trees = tune()
##   min_n = tune()
##   tree_depth = tune()
##   learn_rate = tune()
##   loss_reduction = tune()
##   sample_size = tune()
## 
## Engine-Specific Arguments:
##   lambda = 0
##   nthread = ncores
## 
## Computational engine: xgboost

Note que neste caso nao definimos um grid pois não utilizaremos o tune_grid. O tuning deste modelo será realizado utilizado otimização bayesiana iterativa com a função tune_bayes().

doParallel::registerDoParallel(ncores)

set.seed(321)
xgboost_tune <-
  workflow_adult_xgb_model %>%
  tune_bayes(
    resamples = adult_vfold,
    param_info = xgboost_params,
    # initial = ?,
    iter = 10,  # local=10, kaggle=100
    metrics = metric_set(roc_auc, accuracy, precision, recall),
    control = control_bayes(no_improve = 5,  # local=5, kaggle=30
                            save_pred = T, verbose = F)
  )

doParallel::stopImplicitCluster()

Como cada métrica se saiu nos ajustes:

autoplot(xgboost_tune)

Melhor combinação:

show_best(xgboost_tune, "roc_auc")
## # A tibble: 5 x 13
##    mtry trees min_n tree_depth learn_rate loss_reduction sample_size .iter
##   <int> <int> <int>      <int>      <dbl>          <dbl>       <dbl> <dbl>
## 1     1   143     9         10    5.72e-2 0.000332             0.322     4
## 2     9   254    22         15    1.71e-2 0.000000000145       0.228     2
## 3    16   384    11         13    1.78e-8 0.274                0.208     9
## 4     9   417    25         15    1.87e-7 0.0000408            0.914     8
## 5     8   290    22         10    7.66e-5 0.00000000309        0.863     7
## # … with 5 more variables: .metric <chr>, .estimator <chr>, mean <dbl>,
## #   n <int>, std_err <dbl>
xgbost_best_auc <- select_best(xgboost_tune, "roc_auc")
xgbost_best_auc
## # A tibble: 1 x 7
##    mtry trees min_n tree_depth learn_rate loss_reduction sample_size
##   <int> <int> <int>      <int>      <dbl>          <dbl>       <dbl>
## 1     1   143     9         10     0.0572       0.000332       0.322

Finalizar workflow:

xgboost_workflow_final <- finalize_workflow(
  workflow_adult_xgb_model,
  xgbost_best_auc)
xgboost_workflow_final
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: boost_tree()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 10 Recipe Steps
## 
## ● step_rm()
## ● step_mutate()
## ● step_rm()
## ● step_string2factor()
## ● step_log()
## ● step_normalize()
## ● step_zv()
## ● step_center()
## ● step_scale()
## ● step_dummy()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Boosted Tree Model Specification (classification)
## 
## Main Arguments:
##   mtry = 1
##   trees = 143
##   min_n = 9
##   tree_depth = 10
##   learn_rate = 0.0571636269562865
##   loss_reduction = 0.000332409824875386
##   sample_size = 0.321678699139473
## 
## Engine-Specific Arguments:
##   lambda = 0
##   nthread = ncores
## 
## Computational engine: xgboost

Atributos mais importantes no ajuste do modelo:

xgboost_workflow_final %>%
  fit(training(adult_split)) %>%
  pull_workflow_fit() %>%
  vip(geom = "col")

Finalizar modelo:

xgboost_final <- last_fit(xgboost_workflow_final, adult_split)
collect_metrics(xgboost_final)
## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.860
## 2 roc_auc  binary         0.918

4 Comparar modelos

Como a métrica ara avaliação da competição será a Area Under Receiver Operating Characteristic Curve, vamos avaliar as curvas ROC junto às medias AUC:

bind_rows(
  rpart_final %>%
  collect_predictions() %>% 
  mutate(id = "rpart")
  ,
  xgboost_final %>%
  collect_predictions() %>% 
  mutate(id = "xgboost")
) %>% 
  group_by(id) %>% 
  nest() %>% 
  ungroup() %>% 
  mutate(roc = map(data, ~roc_curve(.x, truth = resposta, `.pred_>50K`)),
         auc = map_dbl(data, ~roc_auc(.x, truth = resposta, `.pred_>50K`) %>% 
                         pull(.estimate) %>% round(4)),
         id = paste0(id, " auc: ", auc)) %>% 
  select(-data) %>% 
  unnest(cols = c(roc)) %>% 
  ggplot() +
  aes(x = 1 - specificity, y = sensitivity, color = id) +
  geom_path() +
  geom_abline(lty = 3)

5 Submissão

Preparar dados para submissão:

adult_val <- readr::read_rds("input/adult_val.rds")

xgboost_model_final <- adult_xgb_model %>% 
    finalize_model(xgbost_best_auc)

adult_fit <- 
  fit(xgboost_model_final,
    formula = resposta ~.,  
    data = bake(prep(adult_recipe), new_data = adult))

adult_val$more_than_50k <- 
  predict(adult_fit, 
          bake(prep(adult_recipe), new_data = adult_val),
          type = "prob")$`.pred_>50K`

Matriz de confusão do pacote caret:

mais informações sobre as métricas na documentação do pacote

adult_val %>% 
  transmute(resposta = factor(resposta, levels = c(">50K", "<=50K")), 
            more_than_50k = ifelse(more_than_50k > 0.5, ">50K", "<=50K") %>% 
              factor(levels = c(">50K", "<=50K"))) %>% 
  table() %>% 
  caret::confusionMatrix()
## Confusion Matrix and Statistics
## 
##         more_than_50k
## resposta  >50K <=50K
##    >50K   2477  1369
##    <=50K   776 11659
##                                          
##                Accuracy : 0.8683         
##                  95% CI : (0.863, 0.8734)
##     No Information Rate : 0.8002         
##     P-Value [Acc > NIR] : < 2.2e-16      
##                                          
##                   Kappa : 0.6144         
##                                          
##  Mcnemar's Test P-Value : < 2.2e-16      
##                                          
##             Sensitivity : 0.7615         
##             Specificity : 0.8949         
##          Pos Pred Value : 0.6440         
##          Neg Pred Value : 0.9376         
##              Prevalence : 0.1998         
##          Detection Rate : 0.1521         
##    Detection Prevalence : 0.2362         
##       Balanced Accuracy : 0.8282         
##                                          
##        'Positive' Class : >50K           
## 
submissao <- adult_val %>% select(id, more_than_50k)
readr::write_csv(submissao, "submissao.csv")

6 Referências

Alguns links úteis na contrução deste kernel: